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We investigate the stability of the pressure-driven, low-Reynolds flow of Brownian sus- 
pensions with spherical particles in microchannels. We find two general families of sta- 
ble/unstable modes: (i) degenerate modes with symmetric and anti-symmetric patterns; 
(ii) single modes that are either symmetric or anti-symmetric. The concentration profiles 
of degenerate modes have strong peaks near the channel walls, while single modes dimin- 
ish there. Once excited, both families would be detectable through high-speed imaging. 
We find that unstable modes occur in concentrated suspensions whose velocity profiles 
are sufficiently flattened near the channel centreline. The patterns of growing unstable 
modes suggest that they are triggered due to Brownian migration of particles between 
the central bulk that moves with an almost constant velocity, and highly- sheared low- 
velocity region near the wall. Modes are amplified because shear-induced diffusion cannot 
efficiently disperse particles from the cavities of the perturbed velocity field. 



1. Introduction 

Microfluidic devices operate in overwhelming laminar conditions where particles mi- 
grate across streamlines either through Brownian motion or shear-induced diffusion 
(SID). In microfiltration, a process to remove unwanted particles from fluid using a mem- 
brane, SID by very s trong crossflow is essen tial to lessen the growth of particle layer cake 
over the membrane ( Vollebregt et al. 2010t ). Particles with differen t sizes are segregated 



in particle fractionation devices due to SID ( Kromkamp et al. 20061) . Cer amic or metallic 



parti cle injection moulding is another subject for SID to play a role in (jKauzlaric et al 



20111) . While Brownian random walk due to thermal fluctuations operates in all times 



the efficiency of shear-induced migration depends on particle-particle and particle-fluid 
interactions. For instance, spherical par ticles move to regions with lower shear rates when 
the P eclet number is sufficiently large (jSemwogerere et al"1l2007t ISemwogerere fc Weeks! 
20081) , but platelike p articles do not sense the details of velocity profile and respond only 



to average shear rate dRusconi fe Stonell2008l) . The lowest reported volume fraction that 
supports SID is < 0.04 dBrown et al.l 1200^1 Well below this limit, SID is turned off as 
the experiments of iRusconi fc Stond (|2008l ) show no transfer of spherical particles from 
the sample to buffer stream of a T-sensor. Interesting and largely unexplored dynamics 
emerges when particle migration is governed not by a single mechanism, but through the 
interplay between Brownian motion and SID. 

Despite the apparent stability of low-Reynolds flows, some trans ient substructures 



like ripples and sharp near- wall features of concentration profiles (see lFrank et al.ll2003 
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Semwogerere et al.ll2007t) are observed in experimental data. Substructures can be long- 



lived or unstable 'modes' excited by anomalous particle migrations, wall roughness, and 
particle-wall interactions. Instabilities in microchannels, however, are very hard to de- 
tect experimentally, and their theoretical prediction is a challenging problem because 
particle migration is generally a slow process compared to the time scale of veloc- 
ity fluctuations. Moreo ver, none of the three kn own instabilitie s induced by inte rfaces 
dHelton fc Yagerll2007lb massive sedimentation (lYiantsiosI [20061; iRao et al.l I2007T) . and 
gravity ([Govindaraian et al.ll2001t ICarpen fc Bradvl 120021 ) seem to occur in microchan- 
nel flows whose streaming velocity profiles are symmetric with respect to the channel 
centreline. Important questions regarding the flow of suspensions in microchannels thus 
include: (i) if excited, how long can stable modes survive, and in what flow conditions 
are they detectable? (ii) what are the shapes of long-lived substructures and how do they 
depend on dimensionless Reynolds and Peclet numbers? (iii) How do Brownian diffusion 
and SID compete in the bulk and near the walls, and when do they collaborate to desta- 
bilize suspension flows? In this paper we attempt to answer these questions, which have 
remarkable implications for the design and manufacturing of microfluidic devices. 

Dynam ics of suspension flows is d escri bed by differen t mod els. The first model intro- 
duced by iLeighton fc Acrivos ( 1987 ) and Phillips et al. ( 1992 ) is phenomenological and 
involves the diffusion fluxes of particles due to particle collisions and the spatial variation 
in the viscosity. The second model roots from the conservation equations of mass, momen- 
tum a nd energy for the fluid and particle phases (|Nott fc Bradvlll994l ). iMorris fc Boulav 
(1999) take into account normal stress diffe rences to hand le curvilinear flows. In this pa- 



per we adopt the constitutive model of lPhillips et al.l (J1992J) for two reasons: (i) combining 
the effects of Brownian and shear-induced diffusions is a straightforward superposition 
(ii) The free diffusion flux parameters allow for an exploration of different flow regimes. 
The literature als o includes models where particle and fluid pha ses interact only through 
Stokes drag (e.g., Rudvak et al. 1997 : Klinkenberg et all 2011 . and references therein). 
However, our study is different in nature: the Reynolds number is smaller than these 
studies by several orders of magnitude, and despite a strong coupling between the parti- 
cle and fluid motions, particles undergo direct two-body collisions governed by shear and 
viscosity gradien ts. These collisions le ad t o particle phase pre s sure a nd shear overlooked 
in the models of iRudvak et al.l (jl997l ) and iKlinkenberg et al.l (|201ll) . 

We assume that the mean streaming velocities of particles and the solvent are identi- 
cal, i.e., the slip velocity is zero because the drag force is high. We include the Brownian 
diffusion in the flux vector. This leads to new steady-state concentrations that do not 
saturate at the centre of the channel. We briefly review the governing equations of sus- 
pension flows in <J2] and find their steady-state solutions in the presence of Brownian 
diffusion. We linearly perturb the diffusion and momentum equations in the vicinity of 
steady state solutions and utilise the Chebyshev tau method to determine the eigenmodes 
of perturbed equations. We present the results of our modal analysis in Sj3] and explain 
the physical mechanism of a new instability that emerges in concentrated suspension 
flows. 



2. Governing diffusion and momentum equations 

We aim at modeling the diffusion inside a microchannel of the width 2W as shown 
in figure a) . We confine our study to regions far from the edges where the flow has 
a two-dimensional nature in the (x, z) plane, which is spanned by the unit vectors e x 
and e z . The x and z axes are along the channel width and flow direction, respectively. 
We define the mean streaming velocity as v = v x e x + v z e z and assume that streaming 
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field remains invariant by changing the y-coordinate. This is a legitimate assumption 
because SID is controlled by shear gradient in the shortest direction. We define cf> as the 
actual concentration of particles, and set its maximum achievable value to <p m = 0.68. 
We scale all lengths by W and all velocities by the maximum velocity V p of the associated 
Poiseuille flow when <fi = 0. Physical quantities are therefore normalized as 4> = <f>/(fi m , 
(x,z) — (x/W, z/W), (v x ,v z ) — (v x /Vp,v z /V p ) and t = V p t/W, where t is the actual 
time and — 1 ^ x ^ +1. From here on, we will drop the bar sign for brevity and will 
explicitly mention if we use actual values. 

For a flow with the mean streaming velocity v, the volume fr action (j) evolves acco rding 
to the following nonlinear partial differential equation (PDE) (jPhillips et al. 19921 ) 

<p,t + (v ■ V)0 = -eV • J, e = cj) m K c (a/W) 2 , 

J= -<f>V (T<P) - p<j?Trr\<&<l> - DV<f>, D = D W/(cj> m K c a 2 V p ), 



(2.1) 



a (with a = 1.82) is the relative viscosity 

(2.2) 



■j 2 

z,x 



2v x 



I 1/2 



where Jis the flux of particles, rj(cj)) = (1 
of the suspension, and 

T = [\4v XiX v z>z - v 2 x , 

is the magnitude of the local shear rate. T is the second invariant of the rate of strain 
tensor. Throughout this paper, (.). s denotes the partial differentiation operator d(.)/ds. 
a is the typical radius of spherical particles in the suspension, and Dq = IcbT / (6irrj s a) 
is the coefficient of Brownian diffusion where rj s and T are, respectively, the solvent vis- 
cosity and temperature, ks is Boltzmann's constant. The coefficients of diffusion fluxes 
K c and K n = j3K c are 'phenomenological' constants. K c and K v indicate the strength of 
two-body interactions due to the spatial variations of collisio n frequency and su spension 
viscosity, respectively. In the absence of Brownian diffusion, Merhi et al. ( 20051) experi- 
mented two flow geometries, parallel-plate and Couette flows, and concluded that j3 > 1 
is independent of flow geometry. Although they use f3 = 2.1 to fit the steady concen- 
tration profiles in both ge ometries, their numer ical results are satisfactory only for the 



Coue tte flow. The model of l Phillips et al.l (1992) performs well for Poiseuille flow studied 
here jStickel fc Powellll2005l section 5). 



One can define the intrinsic time scale is of the suspension flow based on the Brownian 
motion of particles across the channel width described by (a; 2 ) — 4Dt. Since — 1 ^ x ^ 
+ 1, we set the mean square displacement (x 2 ) to 2 2 and define ts = (x 2 )/(4I?) = l/D. 
In H'2.21 ts will be used to quantify the period of transi ent oscillations. For sus pensions, 
the elements of the stress tensor T = [Ty are given as dCarpen fc Bradvl f2002. equation 
2.46) 



Tij = -Sijp + rj{4>) Vv + (V«) 



(2.3) 



where the superscript T denotes transpose. The suspension pressure p is a superposition 
of the solvent and particle-phase pressures. Equation (12.31) is obtained from the stress 
tensor of rheological models (e.g.. lYurkovetskv fc Morrial2008l ) by neglecting differences 
between normal stresses. This assumption is valid when the migration of particles occurs 
in the shear plane, as in the straight channels studied here. 

Defining Re = pWV p /rj s as the channel Reynolds number in the limit of <f) — > 0, the 
normalized continuity and momentum equations read 



V • v = 0, V • T = Re [v.t + (v ■ V) v] . 



(2.4) 



Since Re is very small (10 2 < Re < 1) in most microchannel experiments (ISemwogerere et al 
l2007tlRusconi fc StondE)08r K it is reasonable to work with a constant suspension density 
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Figure 1. (a) Channel geometry; big arrow indicates flow direction, (b) Variation of Pe versus 
<^buik for several choices of Re and /3. We have set K c = 0.03. (c) Steady-state solution of 
the normalized particle c oncentration <f> for /3 = 5. Circles show the experimental results of 
ISemwogerere et al.l (|2007l . figures 3). (d) Solid line: the velocity profile v z o(x) corresponding to 
4>o(x); dashed line: the velocity profile of Poiseuille flow in the limit of <j> — 0. Model parameters 
are K c = 0.03, 0o(±l) = O.21/0 m , W = 25pm, a = 0.7/xm, V p = 2.1 x 10^ 3 m/s, T = 22°, 
r/s = 2.18 x 10~ 3 Pa.s, p = 1232.5kg/m 3 and D = 8.4 x 10" 2 . The actual velocity is computed 
as VpV z o{x). Using these parameters, we obtain Q = 16.92 nl/s, Pe — 137 and 0buik = O.25/0 m . 



p and drop terms like (j)Ap/p where Ap is the difference between the particle and fluid 
phase densities. After solving (|2 .4[) for the velocity field, the particle Peclet number is 



determined as Pe — a 2 V^ z , max / (D W). 

2.1. Steady-state solutions 

In a steady-state fully developed flow, the gradients of physical quantities are nonzero 
only in the x direction, and the maximum velocity occurs at x = (channel centreline). 
The steady solutions <f) = <f) (x), v x = and v z — v z o(x) of equations ([2.1 jl and ()2.4j) are 
obtained by setting V • T = and J = 0. The first relation gives T = 2 |a; | / r^o (^), and 
the latter yields the first-order ordinary differential equation 

^ = -2sign(x)0 2 (l - 4>T [l>+2|a#(l - + 2a(/3 - l)|x|0 2 (l - , (2.5) 

where 770(2:) = , n{4>o(x)). We numerically integrate equation (|2.5p and obtain = <po(x). 
The corresponding profile of v z0 (x) is then determined from v z o(x) = 2£{rio(t;)]~ 1 dl; 
with the boundary condition u z o(±l) = 0. Our numerical experiments show that for large 
values of /3 the concentration 4>q(x) develops a cusp as x — > 0^ and its tails become flat 
for |jc| — > 1. The velocity profile is flattened near the channel centreline by decreasing 
/3. The profile of v z q(x) approaches to 1 — x 2 (v Zymax —> 1) as — > 0. The physical 
values of K c and /3 are determined through fitting the computed profiles of 4>o(x) and 
Vzo(x) to experimental data. In this paper we explore the perturbations of models with 
1 < /3 < 5. In the limit D = the steady solutions are obtained from x(j)o(x)rj l3 ~ 1 = 



constant ((Phillips et al.lll992l ) while </>o(0) is always saturated to 1. We are not interested 
in this extreme case. 

We have used the solution of (j2.5[) and plotted Pe versus </>buik in figure [TJ 6) for sev- 
eral choices of Re and /3. Figures UK c) andQJd) show the steady-state solutions of a 
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suspension flow with j3 = 5. The materi al and geometrical param eters — given in the 
figure caption — of this example come from Semwogerere et all ( 20071 ) for the solvent cy- 
clohexylbromide/decalin mixture. Using an evolution parameter, they have confirmed the 
steady state condition of the flow. We find tftbuik = O-25/0 m and Pe = 137, which agree 
with the experimental data within 5% (see Semwogerere et ai1l2007l figures 3 and 10c). 
T he solid line in figure [If c) accurately reproduces the concentration pro file in figure 10(c) 



of ISemwogerere et al.l (|2007[) whose own analytical predictions (see also Morris fe Boulavl 
1999; iMiller fc Morrisll2006f) show remarkable deviations from the measurements at the 
fully developed stage. This can be due either to electrical stresses, or particle random 
walks. The impressive match between the results of the steady-state model (|2 .5[) and 
experimental data is mainly due to D ^ and supports the latter possibility. The signif- 
icance o f having a non - zero D had not already been discussed/explored in the original 
work of IPhillips et al. (1992) because they worked with very large Peclet numbers of 
0(1O 5 ). 



2.2. Perturbed equations and eigenvalue problem 

The steady concentration of particles in microchannels can easily be disturbed. For in- 
stance, unavoidable surface roughness, sedimentation, particle-wall and particle-particle 
interactions near a wall can induce small amplitude fluctuations on the boundary values 
of </> a- n( i generate global modes. To understand the transient response of suspensions, we 
carry out a linear stability analysis by perturbing the concentration and velocity fields as 
4> = <po(x) + 4>(x, z, t) and v = v z o(x)e z + v(x, z, t). The magnitude of the shear rate and 
the normalized viscosity then become T — Tq(x) +T(x, z, t) and i] = r/o(x) + ^^(^(x)) </>, 
where T = sign.(dv z o/dx) (v x>z +v z , x )- It is remarked that the transient response can 
develop structures in the y-direction as well. Since the steady-state quantities do not de- 
pend on y, the linear response in that direction will include a simple harmonic waveform 
with a wavelength £ y . The magnitudes of all gradients in the y-direction are determined 
by l/£ y . This study is conducted in the long wavelength limit l y —> oo. 

The perturbed equations are simplified by assuming the stream function ip = ipo(x) + 
ip(x,z,t) to express the velocity field: v — (ip, z ,— ip,x)- This implies v (x) — (0,— V>o,x) 
and v = (ip,z, — ip,x)i an d the continuity equation is satisfied automatically both in the 
steady and perturbed states. We now take the curl of (|2.4|) and remove the pressure p 
from computations. The resulting equation together with (|2.1[) are linearized to obtain: 

C n i> + Ci20 = Re (V 2 ^ + A.xxxi - V 2 ^, z ) , (2.6) 
-eV • (J x e x + J z e z ) = 4> >t + (f>o, x ib, z - ipo,x4>,z, (2.7) 

where the perturbed components of the flux vector defined as J x = £21^ + £22</ ) an d 
J z = C311P + Cz2<t>- The linear operators Cij are functions of ipo(x), <fio(x) and their 
x-derivatives. They are obtained by evaluating 

£n = (VxV-T)^, £i 2 = (VxV-T)^ £21 - (J- e x ) ^ , 
C 2 2 = {J-e x )^, £31 = (J- e z ) ^ , £ 32 = (J- e z )j, 

at (ip,4>) — 0. The partial derivatives djdv and d/dg are noncommutative over the 
extended space (v,g) when v = (x,z,t) and g = (^,0) are independent and dependent 
variables, respectively. We have applied the rule [/(t / )5(^)]^g = fd/dv + df jdv to calcu- 
late the partial derivatives in (|2.8|) : i.e., we first differentiate with respect to independent 
variables, then perform the partial differentiations d/dip and d/d<f). The boundary con- 
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ditions associated with the perturbed equations are v(±l,z,t) = and J X (H, z,t) = 0. 
It is remarked that J z (±l) can vary arbitrarily. 

We consider ip(x, z, t) = exp(ifcz — iut)^(x) and </>(x, z, t) = exp(ifcz — \uit)^{x), where 
w = f2 + £i(i = \J — 1) and l z = 2ir/k is the wavelength of oscillations along the channel. 
ft = 2n/(\tB) = 2nD/X is the wave frequency and £ is the growth/decay rate. Short- 
period transient oscillations with A <C 1 are dissolved by thermal fluctuations. Therefore, 
only long-period oscillations of A > 1 can exist. The linear solutions are decoupled in 
the fc-space, and equations (|2.6j) and (|2.7[) remain invariant under the transformation 
x — > —x if ip(x,z y t) — =pip(—x,z,t) and <j){x,z,t) — ±<f>(— x, z, t). This means that both 
symmetric and anti-symmetric modes are supported by the governing equations, and 
that degenerate pairs may exist in the eigenspectrum of uj. Equation (|2.6j) reduces to 
Orr-Sommerfeld stability equation when tj> = 0. 

Most microchannels have typical widths of 2W ~ 10~ 4 m. For a ~ 10~ 6 m, we will 
get e ~ O(10 -5 ) and 0.01 < Re < 1. On the other hand, the response time of <j> to 
particle migrations is scaled by e, which is understood from equation (|2.7|) . Therefore, 
the growth/decay rates of modes supported by diffusion, and not by the perturbations of 
streamlines, can be as small as £ ~ 0(e) ~ 10 -5 , which requires a sophisticated numerical 
procedure to be resolved. 

We utilise the Chebyshev tau algorithm (Orszag 1971 ; Dongarra et al.|[l996f) to com- 



pute the wave functions ^(x) and $(x), and assume 4 f (x) = J2PnT n (x) and $(x) = 
q n T n (x) where T n (x) are Chebyshev polynomials defined over the region —1 < x ^ +1. 
There are six boundary conditions associated with ij) and <fi. This suggests to assume six 
new unknowns, the so-called r variables. Introducing a complex vector z, which contains 
the variables p n , q n (n = 0, 1, • • • , N) and tj (J = 1, 2, • • • , 6), and the Galerkin weighting 
of (|2.6[) and (12.71) . leave us with the linear eigensystem A ■ z = ujB ■ z where A and B 
are complex matrices. In our implementation of the Chebyshev tau method, we use the 



formula d m T n (x)/dx m = 2 m " 1 n(m - l)!C™_ ro (x) (|Gradshtevn fc Rvzhikl 120071 ) where 
C™_ m (x) are Gegenbauer polynomials. We find the generalized complex eigenvalues to 
and their associated right eigenvectors using the routine zggev.f of LAPACK library. We 
begin our calculations with N = 20 and increase N until min[|0|, converges within 
0.5% for the mode with the smallest |£| in the spectrum. The major source of errors is 
the numerical evaluation of the inner products (T n i,CijT n ), especially when the deriva- 
tives of 4>o(x) and To(x) with jump discontinuities at x = appear in the integrand. 
We compute the inner products using a mid-point rule to avoid x = 0, and use finer 
grids in the x-domain as n or n' increase. A uniform grid helps us simultaneously resolve 
the strong near-wall features of certain modes and handle the central cusp. Reaching to 
0.5% error threshold occasionally needs N > 80 to capture short wavelength modes when 
C ~ 0{e). To calibrate our co de, we have s olved the Orr-Sommerfeld stability equation 



and reproduced the results of lOrszagl (|1971l ) up to the 8th decimal place. To assure that 



the physical eigenfrequencies are not sensitive to the choice of basis functions, we used 
Fourier series to reconstruct $(x) and ^(x), and compared the resulting spectra with 
Chebyshev tau algorithm. The results of two methods match very well for modes that 
peak near the channel centreline. The Chebyshev tau method, however, gives more ac- 
curate results for modes that develop bumps near the walls. Our numerical calculations 
show that Fourier series increase the number of spurious modes. 



3. Long-lived and unstable modes 

We first carry out the stability analysis for the steady f3 = 5 model of i j2.ll (see figure 
[TJ. The main properties of this model, which we call model A, are (i) relatively low </>buik! 



Instabilities of Brownian suspensions 



7 



Mode 


D 


Pe 


Re 


</>m</>bulk 


k 


p 


iu n = n + (i 


uj n = fi + Ci 


S1,S2 


0.084 


137 


0.03 


0.24 


1.0 


5 


0.017958 - 0.010418 i 


0.017886 - 0.010543 i 


S3 


0.084 


137 


0.03 


0.24 


1.0 


5 


0.482467 - 0.000397 i 




S4 


0.084 


137 


0.03 


0.24 


1.0 


5 




0.471442 - 0.000425 i 


U1,U2 


0.030 


53.5 


0.08 


0.55 


0.1 


2 


0.000844 + 0.000282 i 


0.000851 + 0.000303 i 


U3 


0.030 


53.5 


0.08 


0.55 


0.1 


2 




0.006675 + 0.000053 i 



Table 1. Eigenfrequencies of long-lived and unstable modes for a suspension flow with 
K c — 0.03 and e = 1.6 x 10~ 5 . Degenerate pairs appear in a single row. The shape of each 
mode is identified through the symmetric ('') or anti-symmetric (J*) shape of &(x). 




Figure 2. The wave functions $(x) (left) and — d$>(x)/dx (right) for the long-lived modes SI 
(top row) and S3 (bottom row). Note the strong near- wall features of mode SI. This property 
is also shared by its symmetric partner S2. Solid and dashed lines correspond to the real and 
imaginary parts of the wave functions. 

(ii) almost no flattening of the velocity profile near the channel centreline. The matrices 
A and B depend explicitly on e and Re. In the limit e = 0, the evolution of (j> is only 
dictated by the deformations of streamlines, and we find only highly-damped, stable 
discrete modes (£ <C —1), which are the characteristics of incompressible Newtonian 
flows at low Reynolds regimes. Turning on the effect of particle migration, e ^ 0, gives 
birth to new long-lived modes with — 1 -C £ < 0. The eigenfrequencies of long-lived 
modes have been calculated for k = 1 and given in Table [1] These modes belong to two 
general families: degenerate and single modes. The oscillation periods and decay rates of 
modes SI and S2 are ps 26 times larger than those of modes S3 and S4. We have plotted 
&(x) and the x-dependent part of v z (x, z, 0) = — exp(ikz)d 1 $(x)/dx in figure [3] for modes 
SI and S3. It is seen that the concentration profile of mode SI has strong peaks near 
the walls while the dominant peaks of mode S3 have been generated near the channel 
centreline. Modes SI and S2 have a better chance for being excited because particle-wall 
interactions can easily disturb the particle concentration and velocity field near the wall. 
Modes S3 and S4 are most likely due to collective random motions because they have 
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Figure 3. The shapes of the unstable modes Ul and U3. Top panels demonstrate the contour 
plots of <j> (left) and v z (right) at t — for mode Ul. Dashed lines in the panels of 4> and v z show, 
respectively, the forms of the associated steady solutions 4>o(x) and v z o(x). The blank region 
corresponds to the lowest 10% of the mode magnitude, which has not been shaded to highlight 
the major near-wall features. Middle and bottom panels show the wave functions &(x) (left) 
and — d^(x)/dx (right) for modes Ul and U3, respectively. Solid and dashed lines correspond 
to the real and imaginary parts of the wave functions. 



small periods of ss 1.1 ts with f2si/^S3 ~ 0(D), decay slowly so that Cs 3 /Csi ~ O(D), 
and have wide-spread patterns. SI and S2 are therefore viscous modes supported by SID. 

We now define the actual decay time of mode X as i% = —(W/V p ) ln(0.1)/£x, which is 
the duration that the amplitude of <j> decays to 10% of its initial value. The actual oscilla- 
tion period of mode X is given by t£ = (W/V p )X x t B - We find (t^if 1 ) = (2.63 s, 4.16 s). 
The characteristic times of mode S3 are quite surprising: (i^j 3 ,ip 3 ) = (69.05 s, 0.15 s), 
which indicate an almost quasi-stationary oscillation in laboratory scales. All these modes 
will exhibit a high signal to noise ratio and can be measured by currently available high- 
speed imagers. The detection of SI and S2 requires high spatial resolution, while S3 and 
S4 need high frame rates due to their small periods. The ripples in the conce ntration pro- 



files, and the near- wall excess /deficit of particles obseryed in experiments (see | Frank et al 



20031: ISemwogerere et al.ll2007t ISemwogerere fc WeekslkoOSl iBrown et alj|2009h might be 



long-lived modes. Since £s3 s» Cs4, we anticipate the coexistence of modes S3 and S4. 
The small difference between their frequencies can yield a quasi-periodic oscillation. 

Increasing 0t>uik has a significant influence on th e velocity profile and flattens it near the 
channel centreline (e.g., Semwogerere et al. 20071 figure 4). Our calculations show that 
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in models with <^>(0) < 1 the flattening of the steady velocity curve is mainly controlled 
by the parameter (3. The flattening of v z o(x) at x = is quantified by the curvature 
C = v z o. xx (0), which equals —2 for Poiseuille flow. To investigate the effect of C on 
the stability of suspension flows, we build another model B, and increase the average 
concentration to </>buik = 0.55/</> m . We then set V p — 6 x 10~ 3 m/s and f3 = 2 to obtain 
a mass flow rate Q = 7nl/s. In this new flow regime, the curvature at x = becomes 
C = —0.02235 indicating a significant flattening. The perturbed equations now result in 
three unstable modes (Ul, U2 and U3) with ( > 0. They have been reported in Table [1] 
for k = 0.1 (£ z = 207r). We have used this particular long wavelength because it expands 
the wavelengths of and ^(x) in the ^-direction, and leads to more visible patterns. 

Unstable modes Ul and U2 are degenerate pairs, and U3 is a single symmetric mode. 
Modes Ul and U2 are the counterparts of SI and S3. The oscillation periods of long 
wavelength unstable modes are: t^ 1 — 31, tp 2 = 30.76 and t^ 3 — 3.92 seconds, and 
their amplitudes are magnified by a factor of 10 within 32-34 seconds for the degenerate 
pair and 181.7 seconds for mode U3. The dominant concentration and velocity peaks of 
modes Ul and U2 are developed near the channel wall (figure[3]). The likelihood of exciting 
these modes is thus very high because of sedimentation or particle-wall interactions. In 
figure[31 we have also demonstrated the contours of Ke[(f>] = Re[exp(ifcz — iut) &(x)] and 
Re[-O z ] = Re[— exp(ifcz — iuit) d 1 $(x)/dx] for mode Ul at t — 0. Prominent curved tails 
have been developed in the wave packets of both <fi and v z . These features are shared 
by Kelvin-Helmholtz type instabilities that are usually triggered at interfaces, but our 
modes have emerged between the central region (where particles move with an almost 
constant velocity) and highly-sheared zone near the walls. 

The channel flow of suspensions with spherical particles involves five parameters (Re, 
Pe, (f>o (0) , f3 and D) and it is impractical to survey the entire parameter space and identify 
its unstable zone. Here we vary only /3, which also controls Pe and the normalized average 
concentration (^buik, and attempt to understand how the gradients of 4>o( x ) and v z q(x) 
near the channel centreline (see §2.1j) correlate with the instability. We find that by 
decreasing f3 in model A, the magnitude of C drops for all modes until mode S3 becomes 
unstable for /3 < 1.8 (Pe « 124.3, bu ik ~ O.29/0 m ). Other modes (SI, S2, and S4) are 
destabilized by further decreasing /3. In another experiment, we gradually increased /? 
in the initially unstable model B to obtain more cuspy concentrations and less flattened 
velocity profiles as |a;| — > 0. We find that mode U3 is stabilized for (3 > 4.5 (Pe « 65.4, 
0buik ~ 0.52/</> m ). Modes Ul and U2 resist stabilization until (3 w 4.8. It is seen that Pe 
should decrease though mildly and </>buik should increase to get instability. Nonetheless, 
both stable and unstable modes exist for 53 < Pe < 124, 0.42 < </>buik < 0.8 and 
1-8 S ft ^ 4.8. The only property shared by all unstable systems is the flattened velocity 
profile, which constitutes a 'necessary' condition for instability. 

Varying the wavelength t z — 2Tr/k has also a notable impact on eigenfrequencies. 
Our calculations show that f2 and |£| increase proportional to k in both the stable and 
unstable models, i.e., unstable modes grow faster for shorter wavelengths. Moreover, we 
find that the wavelengths of $(x) and — d^(x)/dx are approximately proportional to £ z , 
and consequently, the wave packets of degenerate modes become more compact near the 
walls for larger values of k (compare the graphs of — d^(x)/dx in figures [2] and [3]) . All 
these suggest that short wavelength instabilities can rapidly ruin the flow structure in 
the vicinity of the walls. Modes with longer wavelengths and lower oscillation frequencies 
can survive far from the walls and be observed experimentally. In a full three dimensional 
excitation with £ y < oo, we anticipate faster growth and oscillation of unstable modes 
because clumps in the y-direction can probably enhance the migration of particles. 

To this end, we argue that unstable modes Ul and U2 are amplified by the Brownian 
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motion of particles. For a relatively low Pe and due to Brownian diffusion, some particles 
can 'leak' from the region with flat velocity curve and move towards the channel walls. 
Such particles will locally increase the volume fraction and viscosity, while the suspen- 
sion velocity drops. Note that the positive bumps of <j> and v z are out of phase in the 
upper panels of figure [3l Particles initially moving close to the channel walls can also 
penetrate to inner regions through Brownian diffusion. When the overall velocity profile, 
v z o{x) + v z (x,z,t), is flat near the channel centreline, SID cannot efficiently disperse 
particles trapped in the cavities of (0, v). Sustained Brownian migration, back and forth 
between highly- sheared and central regions, thus amplifies the concentration and velocity 
anomalies and triggers a mode. This instability mechanism works for mode U3 as well 
because the shape of v z is sufficiently flat in central regions (bottom panels in figure [3]). 
Modes will be damped when SID becomes dominant due to the steady form of T${x). 
In this condition SID will disperse particles participating in perturbations, and facilitate 
their migration towards the channel centre. This is why the transition from unstable 
to stable modes strongly correlates with the value of /3, which controls the flatness of 
velocity curve. 



4. Conclusions 



We showed that the steady-state solutions of the model of iPhillips et al.l |1992) can 
reproduce recent experimental results of Brownian suspensions with spherical particles. 
We calculated the eigenmodes of the corresponding perturbed equations, and found new 
families of long-lived and unstable modes. Unstable modes that we find occur when the 
fully-developed velocity profile is sufficiently flattened near the channel centreline. The 
fastest growing modes appear as degenerate pairs, and they live inside the highly-sheared 
region near the walls. Since our model has not taken particle-wall interactions into ac- 
count, unstable modes seem to be triggered by the transverse Brownian migrations across 
the streamlines of the fully developed flow. Mode amplification in models with flattened 
velocity profiles is mainly due to inefficient shear-induced diffusion that can, in princi- 
ple, disperse particles and force them to move towards the channel centreline. Rapidly 
growing unstable modes destruct the flow structure near the walls, and they can explain 
the experimentally observed excess and/or deficit of particles near the channel walls. 
The instability mechanism that we suggested operates along the shortest direction of 
the channel. Therefore, dynamics along the neglected y-direction may only affect the 
wavelength and period of developing patterns, and not their general shapes. The three 
dimensional response of Brownian suspensions, especially in channels with 2W ~ H, 
requires further exploration. 

We thank Howard Stone for his insightful comments. A.K. thanks the department of 
Mechanical and Aerospace Engineering at Princeton University for their hospitality and 
generous support. We are indebted to anonymous referees whose criticisms helped us to 
substantially improve the presentation of the paper. 
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